using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQAWCE - Computation of Cauchy principal value.
    /// </summary>
    [Serializable]
    public static class DqawceClass
    {
        /// <summary>
        /// Computation of Cauchy principal value.
        /// </summary>
        /// <param name="f">Double precision function to be integrated.</param>
        /// <param name="a">Lower limit of integration.</param>
        /// <param name="b">Upper limit of integration.</param>
        /// <param name="c">The parameter in the weight function.</param>
        /// <param name="epsabs">The absolute accuracy requested.</param>
        /// <param name="epsrel">The relative accuracy requested.</param>
        /// <param name="abserr">The estimate of the modulus of the absolute error.</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="ier">The requested errorcode.</param>
        /// <returns>The computed Cauchy principal value.</returns>
        public static double Dqawce(HardRealFunction f, double a, double b, double c, double epsabs, double epsrel,
                                    ref double abserr, ref int neval, ref int ier)
        {
            double aa;
            double area;
            double area1;
            double area2;
            double area12;
            double a1;
            double a2;
            double bb;
            double b1;
            double b2;
            double errbnd;
            double errmax;
            double error1 = 0;
            double error2 = 0;
            double erro12;
            double errsum;
            double result;
            double[] alist = new double[QuadPackConstants.LIMIT];
            double[] blist = new double[QuadPackConstants.LIMIT];
            double[] rlist = new double[QuadPackConstants.LIMIT];
            double[] elist = new double[QuadPackConstants.LIMIT];

            int[] iord = new int[QuadPackConstants.LIMIT];
            int iroff1;
            int iroff2;
            int k;
            int krule;
            int last;
            int maxerr;
            int nrmax;
            int nev = 0;
            int limit;

            limit = QuadPackConstants.LIMIT - 1;
            ier = 6;
            neval = 0;
            last = 0;
            alist[0] = a;
            blist[0] = b;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            iord[0] = 0;
            result = 0.0;
            abserr = 0.0;
            if ((c == a) || (c == b) || ((epsabs < 0.0) && (epsrel < 0.0)))
            {
                goto _999;
            }

            //  First approximation to the integral.    
            aa = a;
            bb = b;
            if (a <= b)
            {
                goto _10;
            }
            aa = b;
            bb = a;
            _10:
            ier = 0;
            krule = 1;
            result = Dqc25cClass.Dqc25c(f, aa, bb, c, ref abserr, ref krule, ref neval);
            rlist[0] = result;
            elist[0] = abserr;
            iord[0] = 0;
            alist[0] = a;
            blist[0] = b;

            //  Test on accuracy.   
            errbnd = (((epsabs) > (epsrel * Math.Abs(result))) ? (epsabs) : (epsrel * Math.Abs(result)));
            if (limit == 0)
            {
                ier = 1;
            }
            if ((abserr < (((1.0e-2 * Math.Abs(result)) < (errbnd)) ? (1.0e-2 * Math.Abs(result)) : (errbnd))) ||
                (ier == 1))
            {
                goto _70;
            }

            //  Initialization. 
            alist[0] = aa;
            blist[0] = bb;
            rlist[0] = result;
            errmax = abserr;
            maxerr = 0;
            area = result;
            errsum = abserr;
            nrmax = 0;
            iroff1 = 0;
            iroff2 = 0;

            //  Main loop.  
            for (last = 1; last < limit; last++)
            {
                // Bisect the subinterval with nrmax-th largest error estimate.    
                a1 = alist[maxerr];
                b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                b2 = blist[maxerr];
                if ((c <= b1) && (c > a1))
                {
                    b1 = 0.5 * (c + b2);
                }
                if ((c > b1) && (c < b2))
                {
                    b1 = 0.5 * (a1 + c);
                }
                a2 = b1;
                krule = 2;
                area1 = Dqc25cClass.Dqc25c(f, a1, b1, c, ref error1, ref krule, ref nev);
                neval = neval + nev;
                area2 = Dqc25cClass.Dqc25c(f, a2, b2, c, ref error2, ref krule, ref nev);
                neval = neval + nev;

                //  Improve previous approximations to integral and error and
                // *  test for accuracy.
                // 
                area12 = area1 + area2;
                erro12 = error1 + error2;
                errsum += (erro12 - errmax);
                area += (area12 - rlist[maxerr]);
                if (Math.Abs(rlist[maxerr] - area12) < (1.0e-5 * Math.Abs(area12)) && (erro12 >= 0.99 * errmax) &&
                    (krule == 0))
                {
                    iroff1++;
                }
                if ((last > 10) && (erro12 > errmax) && (krule == 0))
                {
                    iroff2++;
                }
                rlist[maxerr] = area1;
                rlist[last] = area2;
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));
                if (errsum <= errbnd)
                {
                    goto _15;
                }

                //  Test for roundoff error and eventually set error flag. 
                if ((iroff1 >= 6) && (iroff2 > 20))
                {
                    ier = 2;
                }

                //  Set error flag in the case that number of interval bisections
                // *  exceeds limit.
                // 
                if (last == limit)
                {
                    ier = 1;
                }

                //  Set error flag in the case of bad integrand behavior at a point
                // *  of the integration range.
                // 
                if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                    (1.0 + 1.0e3 * Double.Epsilon) * (Math.Abs(a2) + 1.0e3 * Double.MinValue))
                {
                    ier = 3;
                }
                // Append newly created intervals to list. 
                _15:
                if (error2 <= error1)
                {
                    alist[last] = a2;
                    blist[maxerr] = b1;
                    blist[last] = b2;
                    elist[maxerr] = error1;
                    elist[last] = error2;
                }
                else
                {
                    alist[maxerr] = a2;
                    alist[last] = a1;
                    blist[last] = b1;
                    rlist[maxerr] = area2;
                    rlist[last] = area1;
                    elist[maxerr] = error2;
                    elist[last] = error1;
                }

                //  Call subroutine dqsort to maintain descending ordering in the list. 
                DqsortClass.Dqsort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);

                //  Jump out of loop.   
                if ((ier != 0) || (errsum <= errbnd))
                {
                    goto _50;
                }
            }

            //  Compute final result.   
            _50:
            result = 0.0;
            for (k = 0; k <= last; k++)
            {
                result += rlist[k];
            }
            abserr = errsum;
            _70:
            if (aa == b)
            {
                result = -result;
            }
            _999:
            return result;
        }
    }
}